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SECTION  1 
INTRODUCTION 


Surface-to-surface,  surface-to-air  and  air-to-surface  tacti¬ 
cal  missiles  of  interest  to  the  U.  S.  Army  are  being  required 
to  engage  in  increasingly  demanding  scenarios.  The  engagement 
requirements  are  expected  to  get  even  more  exacting  during  the 
next  decade  [1].  These  requirements  can  only  be  met  with  major 
improvements  in  guidance  laws.  The  most  difficult  part  of  tacti¬ 
cal  missile  guidance  is  the  terminal-homing  phase.  Terminal 
homing  generates  missile  commands  to  direct  the  warhead  on  a 
desired  target  or  a  specific  region  of  the  target. 

Proportional  navigation  (pronav)  has  been  used  extensively 
for  terminal  guidance  because  of  its  success  in  conventional 
ground-to-ground  and  ground-to-air  engagements.  In  addition, 
pronav  is  implemented  by  directly  commanding  missile  acceleration 
components  proportional  to  outputs  of  a  gimballed  seeker. 

Several  studies  have  shown  that  pronav  is  incapable  of 
meeting  the  guidance  requirements  in  the  late  1980's  and  beyond. 
The  three  main  reasons  why  pronav  will  not  be  acceptable  for 
future  missiles  are:  (1)  improved  accuracy  requirement  in 
conventional  scenarios,  (2)  more  demanding  future  engagements 
(e.g.,  guidance  of  anti-tactical  ballistic  missiles)  and  (3) 
increased  stress  on  inexpensive  seekers,  gyros  and  accelerometers 
(e.g.,  strapdown  seekers).  Therefore,  there  is  a  need  to  develop 
advanced  estimation,  control  and  signal-processing  techniques 
for  Army  tactical  missiles. 
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1.1  SUMMARY  OF  APPROACH 


The  missile  guidance  problem  is  solved  by  decomposing  it 
into  three  separate  and  simpler  problems  (Figure  1-1): 

1.  Missile  guidance  with  known  missile  and  target 
dynamics, 

2.  Estimation  of  target  dynamics,  and 

3.  Adaptive  autopilot  design. 

The  missile  guidance  solution  is  obtained  using  the  singular 
perturbation  extension  to  the  optimal  control  solution.  Target 
trajectory  and  dynamics  are  estimated  by  fitting  an  Autoregressive 
Moving  Average  (ARMA)  model  to  the  available  measurement  and 
updating  the  model  parameters  recursively  in  time.  Stabilization 
of  the  missile  with  respect  to  random  disturbances  and  tracking 
of  the  nominal  trajectory  provided  by  the  optimal  control  solution 
are  achieved  by  an  autopilot  based  on  an  adaptive  lattice  algorithm 


Figure  1-1  A  Modern  Control-Based  Missile  Guidance  Law 
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1.2  RESULTS 


1.2.1  Advanced  Guidance  Laws  for  3urface-to-Alr  Missiles 

A  singular  perturbation  guidance  law  has  been  developed 
for  medium-range  surface-to-air  missiles  [2].  This  guidance 
law  is  a  significant  extension  of  a  previously  developed  guidance 
law  for  short-range  missiles  [3];  in  medium-range  intercepts, 
the  problem  of  energy  management  should  be  addressed  in  addition 
to  homing  guidance. 

The  mathematical  formulation  has  been  simplified  by  intro¬ 
ducing  separation  of  time  scales.  While  time  constants  for 
medium-range  intercepts  are  significantly  different  from  those 
of  short-range  intercepts,  the  principle  of  time-scale  separation 
is  still  applicable.  The  resulting  simplified  optimal  control 
formulation  requires  solution  to  a  set  of  nonlinear  algebraic 
equations  and  to  an  initial  value  problem,  all  well-suited  for 
on-board  real-time  implementation. 

1.2.2  Target  Trajectory  Estimation 

A  recursive  algorithm  for  estimation  of  ARMA  model  parameters 
from  noisy  samples  has  been  developed.  Application  of  this 
algorithm  to  parameter  estimation  problems  has  exhibited  its 
fast  convergence  and  unbiasedness  in  the  presence  of  noise  [4], 
even  with  short  data  records.  The  algorithm  has  two  versions, 
a  Recursive  Maximum  Likelihood  (RML)  form  and  a  Recursive 
Prediction  Error  (RPE)  form,  both  of  which  posses  a  parallel 
structure  that  makes  them  highly  suitable  for  parallel-processing 
implementation. 

1.2.3  Adaptive  Autopilots 

Lattice-form  algorithms  have  been  developed  for  fast,  recur¬ 
sive  identification  and  control  of  time-varying  systems  [5-7]. 
These  algorithms  have  excellent  numerical  properties  and  a  modular 
structure  that  makes  them  suitable  for  on-board  real-time 


3 


implementation.  Since  the  model  parameters  and  the  control 
input  are  updated  at  each  measurement  point,  recursive  lattice 
algorithms  respond  to  changes  in  system  dynamics  faster  than 
conventional,  nonrecursive  algorithms.  This  advantage  makes 
recursive  lattice  algorithms  a  natural  choice  for  incorporation 
in  adaptive  autopilots.  They  enable  stabilization  of  the  missile 
with  respect  to  random  disturbances  with  short  time  constants 
and  tracking  the  nominal  trajectory  indicated  by  the  optimal 
control  solution. 

1.3  REPORT  ORGANIZATION 

This  report  is  organized  as  follows:  Section  2  summarizes 
the  missile  guidance  problem,  Section  3  develops  advanced  guidance 
laws,  and  Section  4  describes  the  algorithms  for  target  dynamics 
estimation  and  adaptive  autopilot.  Conclusions  are  given  in 
Section  5. 


SECTION  2 


MISSILE  GUIDANCE  PROBLEM  DISCUSSION 


The  missile-target  dynamics  are  highly  nonlinear  partly 
because  the  equations  of  motion  are  best  described  in  an  inercial 
system,  while  aerodynamic /control  forces  and  moments  are  repre¬ 
sented  in  missile  and  target  body  axis  systems.  The  linearization 
of  the  nonlinear  equations  of  motion  is  complicated  by  fast 
changes  in  relative  target  and  missile  velocity  vectors.  There¬ 
fore,  simplified  estimation  and  control  procedures  for  linear 
systems  cannot  be  applied. 

Proportional  navigation  (pronav)  and  other  conventional 
guidance  laws  have  been  developed  using  classical  control  methods, 
based  primarily  on  linear  system  formulations.  It  was  alsr 
necessary  to  divide  the  overall  problem  into  guidance  and  auto¬ 
pilot  design  to  reduce  the  order  of  the  problem.  Since  the 
missile  guidance  problem  is  highly  nonlinear,  time  varying  and 
of  high  order,  classical  methods  are  difficult  to  extend  for 
improved  guidance  laws. 

The  modern  control  theory  formulation  can  treat  nonlinear 
systems  with  multiple  inputs  and  multiple  outputs.  Modern  control 
can  also  handle  trajectory  constraints  (e.g.,  stability  and 
energy  management)  and  terminal  requirements  of  small  miss 
distance. 

In  general,  tactical  missile  guidance  has  three  phases, 
boost,  midcourse  and  terminal.  The  boost  and  midcourse  phase 
involve  energy  management,  threat  avoidance  and  navigation. 

The  terminal  phase  is  described  by  the  seeker  locked  onto  the 
target  (or  a  region  on  the  target).  Guidance  commands  in  the 
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terminal  phase  are  generated  to  intercept  the  target.  The  ter¬ 
minal  phase  is  the  most  demanding  part  of  missile  guidance. 

Our  research  is  focused  on  improving  missile  performance  during 
this  critical  phase. 

2.1  MISSILE  SUBSYSTEMS 

Figure  2-1  shows  important  missile  subsystems  which  impact 
guidance  and  control  laws.  Missile  dynamics  are  driven  by  aero¬ 
dynamic  and  propulsion  system  characteristics,  control  inputs 
and  the  ambient  conditions.  Guidance  commands  must  be  derived 
from  sensor  and  seeker  outputs.  Sensors  measure  missile  inertial 
states,  usually  angular  rates  and  translation  accelerations.  The 
seeker  measures  components  of  target  position  and  velocity  rela¬ 
tive  to  missile  fixed-axis  system. 

The  seeker  derives  information  about  the  target  using  infra¬ 
red,  laser,  radar  or  other  techniques  and  is  often  a  major  part 
of  the  missile  cost.  The  seeker  is  subject  to  jamming  and  decoys. 
Good  signal -processing  techniques  are  often  needed  to  minimize 
the  degradation  of  target  information  due  to  countermeasures. 


Figure  2-1  Major  Missile  Components 
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From  the  guidance  viewpoint ,  seekers  are  either  passive  or 
active.  Passive  seekers  measure  the  two  components  of  line-of- 
sight  rates  or  angular  locations  of  the  target  in  the  missile 
axis  system.  Active  seekers,  in  addition,  provide  measurements 
of  range  or  range  rate.  Passive  seekers  are  highly  desirable 
because  of  lower  cost  and  minimum  vulnerability  to  jamming. 

Each  component  in  the  overall  model  exhibits  strong  non¬ 
linear  behavior.  Missile  aerodynamics  and  propulsion  models  are 
nonlinear,  particularly  because  of  fast  changes  in  speed,  angle- 
of-attack  and  altitude.  Sensors  also  have  large  nonlinear  errors 
like  scale  factors  and  bias.  Seeker  outputs  follow  various 
trigonometric  relationships  causing  additional  nonlinear  behavior. 
Table  2-1  shows  typical  actuator  errors.  These  errors  produce 
additional  nonlinearities  in  missile  dynamics. 


2.2  PAST  APPROACHES  FOR  MISSILE  GUIDANCE 

Many  guidance  laws  have  been  used  for  the  terminal  phase 
in  the  past.  Pursuit  and  proportional  navigation  (pronav)  are 
the  two  most  common  techniques.  More  recently,  pronav  has  been 
used  almost  exclusively  in  advanced  missiles.  Table  2-2  compares 
five  classical  approaches  to  missile  guidance.  Pronav  will 
be  our  baseline  technique  because  it  has  been  most  successful 
in  previous  work. 

In  pronav,  the  commanded  missile  acceleration  is: 


a  =  rt  12- 
amc  n  t 


go 


(2.1) 


where  r  is  the  range,  a  is  the  line-of-sight  rate,  t  is 
the  time  to  go  and  An  is  the  navigation  gain.  In  conventional 
pronav  applications,  t  is  approximated  by 


t 


go 


-r/r  , 


(2.2) 


such  that 
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-  Reduces  controller  Bandwidth 

-  May  cause  instability  in  high  gain 
controllers 


Produces  steady  state  errors  in 
commanded  acceleration 


-  Reduces  controller  bandwidth 
•  May  produce  limit  cycling 


-  Reduces  control  effectiveness  at 
large  commanded  acceleration 


u  a  u£  ♦  random 
deflection 


u  *  control  deflection;  uc  *  commanded  control  deflection 


Table  2-2.  Comparison  of  Classical  Guidance  Laws 


GUIDANCE  LAW 

ADVANTAGES 

DISADVANTAGES 

Comr.and- to-Line-of- 
Sight  Guidance 

No  teminal  seeker  required 

•  Very  inaccurate  against  mowing 
targets  and  with  winds 

•  Data  link  required 

Pursui  t 

•  Noise  Insensitive 

•  Easy  to  use  with  strapdown  seekers 

•  Inaccurate  against  moving  tar- 
gets  and  with  winds 

Proportional 

Accurate  against  constant  velocity 
targets 

•  Inaccurate  against  accelerating 
targets 

•  Stability  is  sensitive  to  noise 

Pursuit  ♦  Pronav 

Between  2  and  3  in  terms  of  accuracy 

•  Between  2  and  3 

Dynamic  Lead 

•  Between  2  and  3  in  terms  of  accuracy 

•  Easy  to  use  with  strapdown  seekers 

•  Between  2  and  3 

•  Stability  problems  if  transition 
to  pronav  occurs  when  signifi¬ 
cant  noise  Is  present 

Navigation  gain  of  between  2  and  4  is  considered  desirable. 

To  achieve  this  navigation  gain,  both  a  and  r  should  be 
used,  o  is  available  from  passive  as  well  as  active  seekers, 
r,  however,  is  measured  only  by  active  seekers.  Conventional 
implementations  of  pronav  with  constant  navigation  ratio,  there¬ 
fore,  require  active  seekers.  The  variation  in  navigation  ratio 
for  stationary  or  constant  speed  targets  is,  nevertheless,  small 
even  with  passive  seekers. 

Characteristics  of  pronav  guidance  laws  are  discussed  in 
Appendix  A.  Appendix  A  indicates  that  errors  occur  in  pronav 
guidance  because  of  the  following: 

1.  Instability  of  pronav  guidance  law  prior  to  impact, 

2.  Target  maneuvers  and  variations  in  target  speed, 

3.  Variations  in  missile  speed, 

4.  Missile  dynamics  and  combined  missile/autopilot 
dynamics,  and 

5.  Sensor  errors  and  seeker  saturation. 

Modern  control  and  estimation  methods  based  on  advanced  guidance 
algorithms  can  overcome  the  problems  indicated  above.  The  next 
section  discusses  the  general  modern  control  theory  formulation. 

2.3  OPTIMAL  CONTROL  SOLUTION 

Three  cases  for  optimal  control  are  considered:  (1)  known 
target  maneuvers  and  missile  states,  (2)  known-target  maneuvers 
but  estimated  missile  states  and  (3)  evasive  intelligent  targets. 
The  resulting  numerical  problems  are  shown  for  each  missile 
guidance  solution. 


2.3.1  Dynamics 


The  dynamics  of  a  missile-target  engagement  may  bo  described 
by  dynamic  equations  of  the  form: 

x  =  f(x,u,aT)  ,  x(0)  =  xQ  , 

0  <  t  <  tf  ,  (2.3) 

where  x  is  the  state  vector  consisting  of:  (1)  components 
of  relative  position,  (2)  components  of  relative  velocity,  (3) 
missile  orientation  angles,  (4)  missile  angular  rates,  and  (5) 
actuator  and  sensor  states,  u  is  the  vector  of  control  surface 
deflections  and  aT  is  the  target  acceleration  vector. 

2.3.2  Optimal  Control  Solution  With  Known  Target  Maneuvers  and 
Missile  States 

The  optimal  control  solution  to  missile  guidance  is  based 
on  minimizing  the  following  function  of  states  and  inputs 

J  =  S(x(tf))  +  j  C(x, u,aT)  dt  .  (2.4) 

The  first  term  in  Eq.  2.4  defines  terminal  requirements  for 
target  intercept.  This  term  decreases  as  the  relative  terminal 
distance  between  the  missile  and  target  decreases.  Requirements 
for  relative  terminal  angular  orientation  for  improved  charge 
detonation  may  also  be  included  in  S(x(tf)).  The  second  term 
specifies  a  preferred  missile  trajectory.  This  term  may  be 
used  to  manage  energy,  to  satisfy  seeker  constraints  or  simply 
to  minimiZ3  flight  time.  Minimization  of  flight  time,  for  exam¬ 
ple,  is  achieved  by  setting  £  to  unity  and  S(x(fj))  =  0. 

The  final  time,  tf,  may  be  constrained  but  is  usually  free 
in  missile  guidance  problems. 
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The  optimization  problem  is  solved  by  using  a  Lagrange 
variable  X  and  a  Hamiltonian  H.  The  equations  leading  to 
the  optimal  solution  are: 

H  =  C(x,u,aT)  +  XT  f(x,u,aT)  , 

i- - («)’  -  -  (#- («f  • 

(iff  -  (if  ♦  (lf>  -  o  ■ 

H(tf)  =  0  .  (2.8) 

Eqs.  2.3  and  2.6  are  coupled.  The  control  is  obtained  from  Eq. 

2.7  and  the  final  time  is  obtained  from  Eq.  2.8.  The  initial 
condition  is  defined  for  states  x  and  the  final  condition  for 
variables  X.  Therefore,  the  computation  of  u  requires  solu¬ 
tion  to  a  two-point  boundary  value  problem  (TPBVP)  represented 
by  Eqs.  2.3  -  2.8.  The  optimal  control  can  be  determined  if 
xQ  and  aT(t)  are  known  and  the  cost  functional  is  defined. 

Therefore,  with  known  initial  conditions  and  target  maneuvers, 
modern  control  requires  a  TPBVP  to  be  solved  for  missile  guidance 
commands.  The  missile  guidance  problem  becomes  more  complex  when 
the  initial  condition  must  be  estimated  from  seeker  outputs  and 
the  target  may  perform  evasive  maneuvers. 

2.3.3  Optimal  Control  with  Noisy  Measurements  and  Estimated 
Initial  State 

Table  2-3  shows  outputs  of  passive  and  active  seekers. 

Clearly  not  all  states  are  measured.  The  measurements  are  also 
noisy. 

A  state  estimator  (e.g.,  Kalman  filter)  is  required  to 
determine  the  current  state  vector  needed  for  optimal  missile 
guidance.  The  extended  Kalman  filter  formulation  may  be  used. 

For  n  states,  n  differential  equations  are  required  for  the 


(2.5) 

X(tf)  *  ||  ,  (2.6) 

(2.7) 
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Table  2-3.  Output  of  Missile  Seekers 


KINO 

OUTPUT 

RELATIONSHIP  TO  MISSILE  STATES 

PASSIVE 

(Infrared,  Optical, 

Two  Une-of-Slght  Angles 

Ratio  of  Lateral  Target  Position  Components 
to  Range  (In  Missile  Axis  Systems) 

Passive  Radar) 

Two  llne-of-SIght  Rates 

Time  Derivatives  of  Above 

ACTIVE 

Range 

Sum  of  Square  of  Relative  Position  Components 

(Active  Radar) 

Range  Rate 

Time  Derivative  of  Range 

state  estimate.  In  addition,  n(n+l)/2  equations  must  be  propa¬ 
gated  for  the  covariance  matrix  and  the  Kalman  gain.  Since  the 
number  of  equations  for  the  covariance  matrix  is  much  larger  than 
the  number  of  equations  for  the  state  estimate,  there  is  a 
significant  payoff  in  simplifying  the  covariance  equations. 

In  linear  systems,  the  error  covariance  matrix  is  indepen¬ 
dent  of  the  measurements.  Therefore,  the  estimation  error  does 
not  depend  upon  the  applied  input.  However,  in  nonlinear  systems 
the  measurements  affect  estimation  accuracy.  The  error  in  state 
estimates  may,  therefore,  be  reduced  by  appropriate  application 
of  inputs.  This  leads  to  a  dual  control  formulation,  where  the 
inputs  provide  guidance  as  well  as  improvement  in  estimation 
accuracy.  The  exact  solutions  to  these  problems  are  quite 
complex  and  simplifications  are  necessary  for  real-time 
implementations. 

2.3.4  Optimal  Control  With  Target  Evasive  Maneuvers 

In  the  previous  two  sections,  we  assumed  that  future  target 
acceleration  time  trace,  aT,  is  known.  If  the  target  is 
capable  of  performing  evasive  maneuvers,  the  future  target 
acceleration  depends  upon  the  missile  flight  path.  At  any  time 
point,  the  target  will  perform  the  most  desirable  maneuver  to 
void  the  missile. 

The  assumption  of  target  evasive  maneuvers  leads  to  a 
differential  game  formulation.  The  target  will  determine  its 


commanded  acceleration  by  maximizing  the  same  or  similar  penalty 
functional  that  the  missile  wants  to  minimize.  The  optimization 
problem  may  be  stated  as 


max  min  J  .  (2.9) 

aT  u 

This  optimization  problem  requires  the  solution  to  Eqs.  2.9  -  2.8 
and  the  following  equation: 


(2.10) 


Note  that  the  resulting  TPBVP  is  even  more  difficult  to  solve 
than  the  one  with  known  target  maneuvers  due  to  the  additional 
constraint  given  by  Eq.  2.10. 


2.3.5  Numerical  Requirements  for  Optimal  Guidance  Laws 

Figure  2-2  is  a  flowchart  for  a  missile  guidance  law  based 
on  modern  control  theory.  Note  that  a  state  estimation  step  is 
required.  In  addition,  the  guidance  law  needs  more  information 
than  with  pronav. 

Table  2-4  shows  the  various  optimal  control  formulations 
for  missile  guidance  and  the  resulting  numerical  procedures 
needed  for  their  solution.  The  numerical  algorithms  are  diffi¬ 
cult  to  implement.  Therefore,  simplifications  which  do  not 
compromise  accuracy  are  needed  for  effective  missile  guidance 
laws. 

Table  2-5  indicates  possible  approaches  for  simplifying 
the  numerical  problems  associated  with  the  optimal  control 
solution.  Because  the  basic  missile-target  dynamics  are  non¬ 
linear,  each  of  these  simplifications  must  be  developed  specifi¬ 
cally  for  the  missile  guidance  law.  A  singular  perturbation 
simplification  technique  has  already  been  developed.  This 
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Figure  2-2.  A  Modern  Control-Based  Missile  Guidance  Law 


Table  2-4.  Optimal  Control  Formulations  for  Missile  Guidance 


FORM*.  AT  ION 

OPTIMAL  CONTROL  PROBLEM 

NUMERICAL  procedures 

Known  Target  Maneuvers  and 
Missile  States 

Optimal  Control  Law  with  Terminal 
Constraints  and  Free  Terminal  Time 

Two-Point  Boundary  Value  Prot’em 
with  Path  Constants 

Noisy  Seeker  Measurements 

Kalman  Filter 

Propagation  of  Differential 
Equations;  Computation  of  Deriva¬ 
tives 

Dual  Control 

Stochastic  Two-Point  Boundary 
Value  Problem 

Differential  Game  (Known  State) 

Difficult  Boundary  Value  Protlem 

Target  Evasive  Maneuvers 

Dual  Differential  Game  (Unknown 
State) 

Difficult  Stochastic  Boundary 
Value  Problem  and  Kalman  Filter 

Table  2-5.  Simplification  of  Optimal  Control 
Numerical  Procedures 


NUMERICAL  PROCEDURES 

SIMPlIF ICAT1DN/AR PRGXIMATION 

resulting  numerical  problem 

Two-Point  Boundary  Value 
Problem  with  Path  Constraints 
(Optimal  Control) 

•  Linearization 

•  Separation  of  time  scales 
(singular  perturbation) 

•  Initial  value  problem 

•  Algebraic  problem 

•  Partially  coupled  solution 

Propagation  of  Differential 
Equations;  Computation  of 
Derivatives 
(Kalman  Filter) 

•  Simplify  Kalman  gain  or  compute 
It  off-line 

•  Time  scale  separation 

•  Reduced  number  of  differen¬ 
tial  equations 

•  Partially  coupled  solution 

Stochastic  Two-Point  Boundary 
Value  Problem 
(Dual  Control) 

•  Define  a  class  of  test  inputs 

•  Convert  to  an  optimal  con¬ 
trol  problem 

Difficult  Boundary  Value 
Protlem 

(Differential  Game) 

•  Parameterize  guidance  law 

•  Reachable  sets 

•  Command  constraint  guidance 

•  Direct  solution  for  missile 
guidance  and  target  evasion 

•  Simplified  optimal  control 
solution 

Difficult  Stochastic  Boundary 
Value  Protlar  end  Kalman 

Filter 

(Dual  Differential  Game) 

•  Define  a  class  of  test  Inputs 

•  Convert  to  an  optimal  con¬ 
trol  solution 

technique  is  described  in  [3].  The  singular  perturbation  guid¬ 
ance  algorithm  for  anti-tactical  ballistic  missiles  ( ATBMs )  is 
shown  in  Section  3.  The  next  section  (Section  2.4)  contains  a 
comparison  of  optimal  and  pronav  guidance  in  one  particular 
ATBM  scenario. 

2.4  COMPARISON  OF  OPTIMAL  AND  PRONAV  GUIDANCE 

We  compare  miss  distances  for  an  anti-tank  missile  and 
for  an  ATBM  due  to  fore-aft  target  acceleration.  The  two 
scenarios  are  summarized  in  Table  2-6.  The  ATBM  has  higher 
speed  and  faster  dynamic  response.  The  pronav  gain  in  each 
case  is  3. 

Miss  in  pronav  occurs  because  of  the  inability  of  the  mis¬ 
sile  to  track  slowing  ballistic  missiles  and  in  the  latter  part 
of  the  flight  due  to  missile  instability.  The  target  tank,  in 
the  first  scenario,  increases  miss  distance  by  stepping  on  the 
brakes. 
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Table  2-6  Description  of  the  Engagement  Scenario 


Figure  2-3  shows  miss  distance  as  a  function  of  tank  decel¬ 
eration  (the  maximum  deceleration  is  expected  to  be  0.5  g).  The 
effect  of  instabilities  is  significantly  amplified  in  the  pre¬ 
sence  of  noise.  Therefore,  the  control  input  is  often  set  to 
zero  in  the  unstable  region.  If  the  control  input  is  set  to 
zero  half-way  through  the  instability,  the  miss  distances  of 
Figure  2-3  will  approximately  double. 

Similar  plots  for  the  ATBM  engagement  are  shown  in  Figure 
2-4.  Most  of  this  miss  distance  results  because  pronav  does 
not  use  the  entire  missile  capabilities.  An  optimal  contro1 
law  gives  miss  distance  as  shown  by  the  solid  line.  The  op 
control  law  is  able  to  achieve  this  improvement  because  it  starts 
applying  acceleration  early  in  the  trajectory  (Figure  2.5). 

The  results  presented  in  this  section  show  that  in  one 
ATBM  engagement  scenario,  the  pronav  guidance  cannot  meet 
performance  requirements,  while  singular  perturbation  does 
adequately  well. 
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Figure  2-4  Comparison  of  Pronav  and  Optimal  Guidance 
Laws  for  ATBM 
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Figure  2-5  Comparison  of  Missile  Acceleration  Profiles 
for  Pronav  Guidance  and  Optimal  Guidance 


SECTION  3 


ADVANCED  GUIDANCE  LAWS  FOR  SURFACE-TO-AIR  MISSILES 


This  section  describes  a  model  for  an  advanced  aircraft 
missile  and  develops  an  optimal  control  formulation  for  improved 
guidance  laws.  A  singular  perturbation  extension  to  the  optimal 
control  solution  is  presented.  The  time  constants  selected  in 
the  singular  perturbation  method  are  based  on  physical  variables 
as  shown  later.  The  resulting  solution  involves  a  set  of  non¬ 
linear  algebraic  equations.  Two  formulations  have  been  used  — 
one  based  on  the  kinematic  state  and  the  other  using  the  energy 
state. 


3.1  MISSILE  MODEL 

With  xm,  vm  and  am,  missile  position,  velocity  and 
aerodynamic  acceleration  vector,  respectively,  the  missile 
dynamic  model  in  the  Cartesian  coordinate  system  is 


*  (x,y ,h)T,  vm  =»  <u,v,w)T  and  am  =  [ax,ay ,az]T.  The 
acceleration  component  along  the  velocity  vector  is  defined, 
because  missiles  in  Army  inventory  do  not  .have  thrust  control 
(the  formulation  can  be  modified,  if  necessary). 


or 


1 

V 


r  •  a 
m  m 


(T-D) 

m 


(3.3) 
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(3.4) 


(T-D)V 

u  ax  +  v  ay  +  w  az - — 

m  is  a  function  of  time  and  its  time  history  is  known.  Thus, 

we  can  consider  a.  a  and  a  as  our  control  variables  with 

x  y  z 

the  constraint  of  (3.4). 

To  simplify  the  derivation  of  the  guidance  law,  we  convert 
the  velocity  equations  into  a  total  speed  and  two  flight  angle 
equations.  Equations  3.1  and  3.2,  thus,  become 


where 


=  V  cos  y  cos  4>  , 

=  V  cos  y  sin  <J>  , 

■  V  sin  y  » 

*  (T-D)/m  -  g  sin  y  » 

=  L  sin  a/m  V  cos  y  > 


Y  -  (L  cos  a  -  mg  cos  Y)/m  V 


(3.5) 

(3.6) 

(3.7) 

(3.8) 

(3.9) 


(3.10) 


4>  and  y  are  flight  path  angles  in  the  horizontal  and  the 
vertical  planes  (see  Figure  3-1).  L  sin  a  and  L  cos  a  are 
lift  components  in  the  x-y  and  vertical  planes,  respectively. 
Thus,  the  two  control  variables  are  given  explicitly  in  this 
formulation . 

The  thrust  time  history  is  predefined.  The  drag  is  written 
as  follows 


D  =  |  P(h)  V2s  CA  , 


(3.11) 


CA  -  CA  +  CA  |a|  +  CA  a  , 
o  a 


q  -  \  p(h)  V2  . 


(3.12) 


(3.13) 


a  is  the  total  angle-of-attack  and  the  density  p  is  a  function 


T 


X 


Figure  3-1  Axis  System  for  Missile  Dynamics 

of  altitude.  The  total  angle-of-attack  is  related  to  total  lift 
as  follows 

L  =  |  p(h)  V2s  CN  .  (3.14) 

The  target  position  velocity  and  acceleration  vectors  are  related 
to  each  other 

xT  =  vT  (3.15) 

vT  =  aT  (3.16) 

The  final  condition  for  missile  intercept  is 

xm(tf)  =  xT(tf)  (3.17) 

where  t^  is  free.  The  optimal  control  solution  will  be  based 
on  the  minimization  of 
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f*  f 

=  /  C  ( X.  v, 

Jo  m  m 


am>  dt 


(3.18) 


Typically  £  is  a  function  of  drag,  total  speed,  energy  loss 
rate  and  total  flight  time.  If  the  optimization  criterion  were 
minimum  time,  the  performance  index  would  be 


£  =  1 


(3.19) 


and  for  the  minimum  energy  loss  solution 


£  -  D  •  V 


(3.20) 


In  our  discussion,  it  is  assumed  that  xT(t)  is  known, 
based  on  an  estimate  of  the  target  motions. 


3.2  OPTIMAL  SOLUTION 

The  solution  to  the  optimal  control  problem  is  obtained 
by  defining  a  Hamiltonian  as  follows.  The  Lagrange  multiplier 
method  for  the  optimal  solution  gives  the  following  Hamiltonian 
for  this  problem 

H  =  Xx  v  cos  Y  cos  <f>  +  Xy  V  cos  y  sin  <f>  +  Ah  V  sin  y 
+  Ay[(T-D)/m  -  g  sin  y]  +  ^  L  sin  o/(mV  cos  y) 


+  X^(L  cos  0  -  mg  cos  y)/(mV)  +  £  . 


(3.21) 


Ax,  Ay,  Ah,  Ay,  A^,  Ay  are  the  various  Lagrange  parameters 
and  £  is  a  general  cost  functional.  Since  the  Hamiltonian 
does  not  depend  on  x  or  y 


Ax  =  Ay  =  0  . 


(3.22) 
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The  other  four  Lagrange  variables  are 


xh  * 

-3H/3h  , 

xh(tf ) 

free  , 

(3.23) 

*V  = 

-3H/3V  , 

*v(tf ) 

=  o  , 

(3.24) 

II 

-e- 

•  r< 

-3H/a<}>  , 

W 

=  0  , 

(3.25) 

• 

X  = 
Y 

-3H/aY  , 

W 

=  0  . 

(3.26) 

The  optimality  conditions  give 


9H 

3L 


V  3D 

■  ~m  at  +  x<j>  sin  o/(mV  COS  y) 
+  X^  cos  a/ (mV)  +  «  0  , 


(3.27'' 


3  H 

3^  =  x4>  L  cos  c/(mV  cos  y)  -  Xy  L  sin  a/( mV)  =  0.  (3.28) 
Note  that  neither  D  nor  £  should  depend  on  a.  Thus, 


tan  a 


X  COS  Y 
Y  T 


(3.28) 


Since  the  final  time  is  free, 


H(tf )  *  0  •  (3.30) 

Note  that  H  is  an  explicit  function  of  time  through  T  and 
m.  Thus,  H  is  not  zero  throughout. 

The  unknowns  are  Xx,  Xy,  L  and  o  and  the  time  histories 

of  xh»  XV’  x$»  and  Xy*  Thus,  six  forward  and  four  backward 
differential  equations  need  to  be  solved  in  a  two-point  boundary 
value  problem  (TPBVP). 
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3 . 3  TIME-SCALE  SEPARATION 


The  six  missile  dynamics  equations  are  scaled  as  follows 


r  d(x/r ) 
V  dt 


=  COS  Y  COS 


<t>  » 


(3.31) 


r  d(y/r ) 
V  dt 


=  COS  Y  sin  $  » 


(3.32) 


*\nax 

V 


d(h/h 

dt 


max' 


sin  y  • 


(3.33) 


The  time  constant  associated  with  the  total  speed  equation  is 
determined  by  substituting  for  a  in  terms  of  lift  in  the  drag 
equations  (Eqs.  3.11  -  3.14) 


dV  =  T 
dt  m 


1 

m 


QsCa  + 


+ 


C.  L 

V 


QsC. 


N 


a 


(3.34) 


The  total  speed  equation  has  two  parts.  The  first  term  is  known 
and  can  be  large  when  the  thrust  is  on.  The  second  term  depends 
on  flight  condition  and  control  input  and,  thus,  controls  the 
dynamics  of  the  total  speed  equation.  The  drag  term  varies  from 
fractions  of  a  g  to  a  few  g's  in  most  medium-range  missiles.  The 
gravity  term  is  always  less  than  one  g  (because  w/V  1).  Thus, 
dividing  by  g  the  right-hand  side  becomes  of  unit  order.  The 
time  constant  associated  with  this  equation  could  change  if  the 
missile  drag  was  larger  by  a  factor  at  high  angle-of-attack. 


V  /I  dV\  (T-D)  w 
g  \V  dt/  “  mg  "  V  * 


(3.35) 


Similarly,  a  time  constant  can  be  developed  for  the  total 
energy  equation 
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dE  (T-D)  v 
dt  mg 


(3.36) 


Since  D  is  of  order  of  mg,  the  nondimensional  equation  for 
E  is 


V  1  dE  =  2  (T-D) 

g  E  dt  ,  .  2gh  mg 

V2 


(3.37) 


Thus,  the  total  energy  equation  has  essentially  the  same  time 
constant  as  the  speed  equation. 

The  time  constants  associated  with  $  and  y  are  obtained 
as  follows 


mV  d$  =  _ L_ 

L  dt  L 
max  max 


sin  a 
cos  y 


(3.38) 


mV  dy  _  L 

L  dt  L 
max  max 


cos  a 


+  mg  cos  y 
max 


(3.39) 


Thus,  the  time  constant  associated  with  flight  path  angles  is 


V/%ax  '  <3-40> 

The  missile  lateral  acceleration  commands  are  generated  by  the 
autopilot.  The  particular  separation  of  time  scales  for  short- 
range  missiles  was  discussed  previously.  Table  3-1  shows  that 
the  time  scales  for  missiles  with  significant  cruise  phase 
(medium-range  intercepts)  could  be  significantly  different. 

A  clear  separation  of  time  scales  associated  with  various  states 
is  clear. 

Slowest :  position  components,  x  and  y,  and  total 

speed  equation, 

Slow:  altitude,  h, 

Fast :  flight  path  angles,  and 

autopilot  states. 


Very  fast: 
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Table  3-1  Time  Constants  Associated  With  Various  States 
During  Midcourse 


COMPONENTS 

TIME  CONSTANT 

VALUE  (s) 

Position 

X 

r/V 

40 

y 

r/V 

40 

h 

*wv 

6 

Velocity 

Total  speed  (or  energy) 

V/g 

50 

Flight  path  angles 

V/IW 

max 

1.67 

Accelera- 

L 

Autopilot 

0.6  to 

tion 

Orientation 

time  constants 

0.05 

r  *  20,000  m,  h  *  3,000  m,  V  *  500  msec”*,  a  *=  30  g's 

m  max 


CONTROL  SJRCAC£  DEFLECTIONS 


Figure  3-2  Schematic  of  Singular  Perturbation  Guidance 
Logic  for  Midcourse  Phase 


Note  that  the  distribution  of  time  scales  is  different  from 
that  realized  for  the  end-game  problem. 

The  approach  to  obtain  guidance  commands  using  singular 
perturbation  methodology  during  midcourse  phase  is  shown  in 
Figure  3-2. 

3.4  SIMPLIFIED  OPTIMAL  CONTROL  SOLUTION  BASED  ON  SEPARATION 
OF  TIME  SCALES 

We  will  solve  the  guidance  problem  in  four  parts,  each 
part  corresponding  to  one  of  the  time  scales. 

3.4.1  Slowest  Time  Scale 

x,  y  and  V  are  the  slowest  variables.  Since  h,  4> 
and  y  are  faster,  those  equations  may  be  considered  to  be  in 
equilibrium.  Therefore  (the  variables  in  the  slowest  time  scales 
are  denoted  by  superscript  "1"), 

Y1  =  a1  =  0  ,  (3.41) 

L1  =  mg  .  (3.42) 

Thus,  the  Hamiltonian  simplifies  to 

H1  =  A*  V1  cos  t1  +  A*  V1  sin  +  xJ(T-D)1/m  +  C.(3.43) 
x  y  v 

The  adjoint  equations  for  A^  and  A^  give 

A^-  1 

3H  V  3D1  ,  3C 

3h  u  “  "  m  7h“  3h 

+  (A*  cos  (fr1  +  A l  sin  4>1)(-g/V1)  (3.44) 

x  y 

||  *  0  *  -A*  V1  sin  t1  +  xj  y1  cos  4*1  .  (3.45) 
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because 

V1  =  V2  -  2gh 
Eq.  3.45  gives 

tan  (J)1  =  xj/xj  (3.46) 

y  x 

Since  XX  and  X*  are  known  to  be  constants,  the  flight  path 
x  y 

in  the  horizontal  plane  for  the  slowest  dynamics  is  a  straight 
line.  The  straight  line  must  joint  the  horizontal  projections 
of  the  desired  current  and  final  conditions.  In  addition,  since 
H1(tf)  and  Ay(tf)  are  zero, 

Xx  yl(tf)  cos  4*1  +  V1(tf)  sin  01  +  C1  (t.f)  -  0  .(3.47) 

Eqs.  3.46  and  3.47  give 


X 


1 

x 


-C1(tf)  cos  (jj1 
V1(tf) 


X 


1 

y 


C(tf)  sin  4*1 
V1(tf) 


(3.48) 


(3.49) 


The  adjoint  equation  for  Xy  is 


:i 


X^  i 

AV  3D1 


1C 


Xy  =  K  cos  <T  +  Xj  sin  *  -  +  fv 


3_ 

av 


c1(tf) 

vx(tf) 


x1  1 

3D_  1 

m  av  *  xv^f; 


o  . 


(3.50) 


Eqs.  3.44  and  3.50  give  the  optimal  altitude  in  slowest  time-scale 
approximation.  This  still  requires  the  solution  to  a  first-order 
TPBVP.  This  requirement  could  be  removed  if  the  Hamiltonian  was 
not  an  explicit  function  of  time. 
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3.4.2  Altitude  Dynamics 


To  solve  for  the  optimization  problem  in  the  altitude 

dynamics  time  scale,  Xx,  Xy  and  Xy  can  be  used  from  the 

previous  solution.  Since  flight  path  angle  dynamics  are  still 

•  • 

faster  than  altitude  dynamics,  setting  <fi  and  y  to  zero  we 
get  (variables  in  sltitude  dynamics  time  scales  are  written 
with  supercript  "2") 


a2  *  0  ,  L2  -  mg  cos  y2  . 


(3.51) 


The  Hamiltonian  can  thus  be  modified  to 

H2  *  xj  V1  cos  y2  cos  (j)1  +  V1  cos  y2  sin  Q1 
x  y 

+  X2  V1  sin  Y2  +  Xy[ (T-D2 ) /m-g  sin  y2]  +  C  .  (3.52) 

The  adjoint  equation  X^  gives 

-  0  -  -Xx  V1  sin  y2  cos  01  -  xj  V1  sin  Y2  sin  <j)1 


+  (X2  V1  -  g  cos  Y2  =  0  • 


(3.53) 


C  is  not  likely  to  be  a  function  of  y.  Thus, 
2  (X2  V1  -  g  Xy) 

tanif  °  -  c(f-)  • 


(3.54) 


The  adjoint  equation  for  X^  is 


^h  XV  3D2  3 C  \2(t  )  is  free 

^t  ~  3K  ’  Xh(tf}  is  free  * 


(3.55) 


2  2 
Xh(tf)  may  be  determined  from  the  relationship  that  H  (tf) 

is  zero.  Thus,  the  optimal  flight  path  angle  may  be  determined 

from  (3.54)  and  (3.55). 
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The  approach  to  compute  flight  path  angle  dynamics  has 
already  been  solved.  Reference  [1]  provides  an  approach  to 
direct  the  missile  from  an  initial  flight  path  to  the  desired 
flight  path  angles  (<J>,y).  As  has  been  shown  previously,  lift 
is  applied  perpendicular  to  the  plane  containing  the  initial 
and  the  desired  velocity  vector.  The  magnitude  of  the  lift 
is  proportional  to  the  square  root  of  the  angle  through  which 
the  flight  path  must  be  changed  (see  below). 


3.5  ALTERNATE  ENERGY  STATE  FORMULATION 

The  optimality  equations  and  the  corresponding  simplifica¬ 
tions  have  been  obtained  using  the  energy  state  to  replace  the 
total  speed  state  in  the  missile  equations  of  motion.  This 
offers  certain  computational  advantages  as  will  be  seen  in  the 
optimality  equations.  With  E  replacing  V,  the  state 


equations  are 

x  *  V  cos  y  C°S  <f>  (3.56) 
y  *  V  cos  y  sin  <t>  (3.57) 
h  =  V  sin  y  (3.58) 
E  =  V(T-D) /mg  (3.59) 
<J>  =  ( L  sin  o)/mV  cos  y)  (3.60) 
Y  =  (L  cos  a  -  wg  cos  y)/mV  (3.61) 


The  performance  index  will  be  written  as  a  combination  of  the 
flight  time  and  the  energy  loss  during  flight. 


J 


( 1-0  dt  +  ? 


E(tf  ) 


(3.62) 


Eq  is  a  normalizing  factor  and  represents  nominal  energy  loss 
per  unit  time.  General  performance  indices  can  also  be  studied. 
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The  Hamiltonian  for  this  problem  is 


H  “  Xx  V  cos  Y  cos  <f>  +  xy  V  cos  y  sin  <f>  +  Xfa  V  sin 
+  V(T-D)/mg  +  (X^  L  sin  o)/(mV  cos  y> 

+  Xy(L  cos  a  -  mg  cos  y'i/m'V  +  (1  -  5) 

The  optimality  equations  are 


Xx  -  Xy  *  0 


Xh  -  -  3H/3h 


Xh(tf)  *  free 


X£  «  -  3H/3E 


XE^f^  “ 


X^  =  -  3H/3<j> 


X^Ctj)  =  free 


Xy  »  -  3H/3y 


X  (t{)  =  free 


The  optimality  equations  are 


3H 

3L 


XEV  3D 


X.  sin  o 

A. 


K  d  LI  ,  <p  -  XY  C°S  ° 

mg  3L  mV  cos  y  +  mV 


=  0 


or 


3H 

3a 


X.  L  cos  a  XL  sin  a 
- I - - -  =  0 


mV  cos  y 


mV 


tan  a  = 


Xy  cos  Y 


Since  the  final  time  is  free 


H(tf)  »  0 

Again  a  four  time  scale  solution  is  sought.  The  ordering 


(3.63) 

(3.64) 

(3.65) 

(3.66) 

(3.67) 

(3.68) 

(3.69) 

(3.70) 

(3.71) 

(3.72) 
of  the 
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time  scales  is  as  follows: 


x,  y,  E 
h 

<J>>  Y 

Autopilot 

Dynamics 


slowest 

slow 

fast 

fastest 


3.5.1  Slowest  Time  Scale 


-  °i  -  0 


L,  =  mg 


and 


(3.73) 


Hx  =  Ax  V^cos  +  Ay^V^sin  <J>X  +  A£  V1(T-D)/mg 

+  (1  -  C)  (3.74) 


A  and  A  are  constants.  Optimal  <j>  and  h  are  obtained 
1  yl 

from 

Xy 

=  -  Ax  V1sin  <J)1  +  Xy  V^cos  $x  =  0  =>  tan  «j>1  =  y^-  (3.75) 


£H 

3h 


[A  cos  <j>x  +  Ay^sin  ^1  +  XE1^T_D^/mg^~g/Vl-* 


XE  V 
E1  3D 

mg  3h 


Note  that  since  Vx  =  V 2g(E-h) ,  3Vx/3h  =  -g/Vx 


(3.76) 


(3.77) 


H(tf)  -  0 


gives 
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(3.78) 


(A^cos  ♦1  +  Ay^sin  <J>1)V1(tf)  +  -^mg  [T(tf )-D(tf )] 


+  (1-5)  -  0 


or 


5[T(tf)-D(tf>]  n 

(A  cos  +  A  sin  )  - - - — - - (3.79) 

xl  1  yl  1  ^O1118  V^(tf) 


Eqs.  (3.75)  and  (3.79)  give  the  solutions  for  A  and  A 

X1  yl 

Nominally,  a  differential  equation  must  be  solved  to  obtain 
£.  Because  the  Hamiltonian  is  explicitly  time  dependent  through 
T  and  m  the  Hamiltonian  is  not  constant  throughout.  We  will 
use  an  approximation  which  will  avoid  the  need  to  solve  the 
differential  equation  in  the  backward  direction.  The  approxi¬ 
mation  consists  of  using  average  values  of  T  and  m  in  the 
definition  of  the  Hamiltonian.  Thus, 


Vi^av'0!) 

Hx  =  A^y^os  4>1  +  A^VjSin  ^ — -  +  (1-0(3.80) 


mav8 


With  this  approximation 
H1(t)  =  0 


Therefore, 


X  ^av-V  ,  1-5 

E1  mav8  V 


C[T(tf)-D1(tf)]  1-c 

TT  S  /  i 


mav8E0 


V1(tf ) 


(3.81) 


Using  Eq.  (3.80),  Eq.  (3.76)  becomes 


XE  V1 

H  (1-5)  “  m  If  =  0 

V?  mavg 


(3.82) 
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£[T(tf )-D(tf)  1_£ 


mavgE0 


V1(tf) 


3D 

1  3h(3.83) 


This  equation  may  be  solved  to  obtain  the  optimal  altitude.  Note 
that  the  three  adjoint  variables  are 


cos  <j>^ 


£  [T(t^  )-D^(t^ )  ]  ^ ^ 


E0mg 


Vx(tf ) 


(3.84) 


sin  4^ 


K [T(t j )-D^(tj) ] 

■n  ^  tt  /  i  v 


Eomg 


V1(tf ) 


(3.85) 


m  g 
av& 


^av'V 


^[T(tf  )-D1(tf )]  (l_o  (l-£) 

"  +  ”  '  •  '  -  v 

V1 


E0mg 


V1(tf) 


(3.86) 


The  two  special  cases  corresponding  to  minimum  time  (E,  =  0)  and 
minimum  energy  loss  (£  =1)  are  shown  in  Table  3-2.  For  the 
minimum  energy  loss  case,  the  drag  is  the  lowest  at  the  optimal 
altitude. 

The  flight  path  direction  in  the  horizontal  plane  is  deter¬ 
mined  by  the  horizontal  projection  of  the  line  which  joins  the 
current  missile  position  to  the  intended  final  missile  position. 


3.5.2  Altitude  Dynamics 

To  solve  for  the  altitude  dynamics  ,  X  ,  and  A„ 

X1  yl  E1 

are  used  from  the  previous  time  frame.  Since  flight  path  angle 
dynamics  are  still  faster  than  altitude  dynamics,  setting 
and  y  to  zero,  we  get  the  following  equilibrium  conditions 

a2  =  0  (3.87) 

Lg  *  mg  cos  y2  (3.88) 
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Table  3-2  Optimal  Solutions  in  the  Slow  Dynamics  for  the  Minimum  Time 
and  the  Minimum  Energy  Loss  Problems 


where  the  subscript  '2'  denotes  variables  in  the  altitude 
dynamics.  The  corresponding  Hamiltonian  can  be  written  as 


H2  =  Xx  V2cos  Yacos  4>1  +  Xy  V2cos  Y2sin  <j> 1  +  xn  V2sin  y2 
112 


+  XE^V2(T-D2)/mg  +  (1-5) 


(3.89) 


The  optimal  flight  path  angle  in  the  vertical  plane  is 
3H2 

=  o  =  -(Xx  cos  4> ^  +  Xy  sin  4>1)V2sin  y2  +  Xh  V2cos  Y2 
211  2  (3.90) 

Using  Equation  (3.26),  we  get 


tan  Yn  = 


2  (X  cos  4>.  +  X  sin  (J>-  ) 


(3.91) 


Normally  we  would  have  to  solve  a  differential  equation  to  obtain 
Xh  because  T  and  m  are  explicit  time  functions.  However,  we 
would  again  set  T  and  m  to  their  average  values  which  gives 
an  additional  equation: 


H2(t)  =  0 


(Xv  cos  <|)1  +  Xtt  sin  (ji^VgCos  y0  +  Xh  V9sin  y 


xi  '  1  *1 


2  h_  2‘ 


<Tav“D2) 
+  X  V  — — — — 

E1  2  mav2 


+  (1-5)  =  0 


(3.92) 


Using  Eq.  (3.91)  to  eliminate  X.  ,  we  get 

n2 

XEl(Tav-D2> 

(X  cos  4>-  +  x  Sin  4>-  )sec  y2  +  - = — z -  +  -w- 

X1  1  yl  1  ^  mavg  v2 


=  0  (3.93) 
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Eq.  (3.80)  is  rewritten  as 


(A  cos  4>-,  +  A  sin  <!>.,)  + 


XEl(Tav-Dl> 


[-  i  y-  i  m  g 

1  *  1  av6 


+  i^-0 

V1 


(3.94) 


Subtracting  Eq.  (3.94)  from  Eq.  (3.93) 


XE  <D1'V 


1 

(A  cos  +  A  sin  <j)1)(sec  Y0-l)  +  - - — - 

xi  1  1  2 


(1-5)(V1-V9) 

+  - VT . .  =  0 

vlv2 


(3.95) 


From  Eq.  (3.76)  and  Eq.  (3.80)  we  can  find  one  solution  to  A 


\  (H)g/V^ 

HT^i  *  3D1/3h 


(3.96) 


Therefore  Eq.  (3.95)  is  simplified  to 


(Ax  cos  4^  +  A  sin  (J^Ksec  Y2  - 


V1~V2  (D1-D2>« 

_  (l-E)  — - — —  +  — = — - - 

v  s  '  VV  ? 

V1V2  V~  SDj^/ah 


(3.97) 


This  equation  gives  a  value  of  sec  Yg.  Y2  is  positive  if 

h2  <  ^1  anc*  *s  negative  if  <  hg.  Note  that  if  this  equation 

solves  out  to  a  negative  value  of  (sec  Yg-l),  the  desired 

flight  path  angle  is  ±90°  (this  happens  because  Eq.  (3.90)  must 

be  modified  when  optimal  Yg  is  90°).  Note  that 

(A  cos  +  A  sin  <$-1)  is  obtained  from  Eq.  (3.84)  and  (3.85) 

X1  1  yl 

'5[T(tf)-D1(tf)]  (1  -  ‘ 

V°S  »i  +  xyisln  *1  ’  -[ - E^g-  —  +  tJtvJ  <3'98> 
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3.5.3  Attitude  Dynamics 

Given  desired  values  of  the  flight  path  angles  and  the 
current  values  of  flight  path  angles,  the  lift  vector  and  its 
orientation  can  be  computed  using  techniques  developed  previously 
[8].  We  shall  summarize  them  here  since  they  form  an  important 
part  of  the  overall  guidance  law.  The  orientations  of  the 
current  velocity  vector  and  the  desired  velocity  vector  are 

v  =  [cos  Y  cos  <|>,  cos  Y  sin  <J> ,  sin  y]  (3.99) 

c 

vd  =  [cos  Y2  cos  cos  Y2  sin  $1#  sin  Y2]  (3.100) 

The  total  angle  through  which  the  flight  path  must  be  changed  is 
given  by 

=  arcos  (v *vd)  (3.101) 

The  lift  vector  net  of  gravity  must  be  perpendicular  to  v  , 
in  the  plane  containing  vc  and  v^.  The  value  of  the  lift 
depends  on  A4>  and  the  variation  of  drag  with  lift.  It  was 
shown  previously  that  the  lift  is  of  the  form 

L  %  K  /“Eijr  (3.102) 

K  may  depend  on  dynamic  pressure  and  time-to-go.  The  value 
of  K  may  be  derived  from  the  results  of  [8]. 

3.6  ALGORITHM 

The  algorithm  involves  the  following  steps: 

1.  Based  on  current  missile  state,  compute  y,  $, 
density,  drag  and  total  energy. 

2.  Computer  desired  <[> ,  based  on  the  line  joining 
the  current  missile  location  in  the  horizontal  and 
the  desired  missile  location  at  the  end  of  the 
midcourse  phase. 
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3.  Based  on  the  drag,  compute  an  approximate  value  of 
Vj^Ctj),  D(tf)  and  time-to-go.  Time-to-go  is 

used  to  determine  Tav  and  mav. 

4.  Solve  Eq.  (3.83)  for  the  optimal  altitude.  Eqs. 
(3.84),  (3.85)  and  (3.86)  give  X„  ,  X,  and  X_,  . 

X1  yl  E1 

5.  Solve  Eq.  (3.97)  for  Y2«  Note  that  for  computa¬ 
tional  advantages  (sec  Y2-l)  should  be  approxi¬ 
mated  as 

2  sin2(Y2/2) 
cos  y2 

particularly  when  Y2  is  small. 

6.  Compute  vc,  vd  and  Aip  using  Eqs.  (3.99)  to 

(3.101).  A  numerically  desirable  way  to  compute 
Aip  uses  the  formula 

Aip  =  2  arcsin 

The  vertical  bars  represent  the  2-norm  of  a  matrix. 
Compute  the  desired  lift  using  Eq.  (3.102). 

7.  The  orientation  of  the  lift  vector  is  given  by 
Vc  and  Vd. 

8.  Add  the  component  of  gravity  perpendicular  to  the 
velocity  vector  to  the  lift  computed  in  steps 
(3.59)  and  (3.60). 
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SECTION  4 


TARGET  DYNAMICS  ESTIMATION  AND  ADAPTIVE  TRAJECTORY  TRACKING 


4.1  EQUATIONS  OF  MOTION 

For  state  estimation,  we  will  describe  the  relative  dynamics 
of  the  target  with  respect  to  the  missile  in  the  inertial  axis 
system.  If  x  is  the  vector  of  relative  position,  vt  and  vm 
are  vectors  of  target  and  missile  velocity,  respectively,  and 
at  and  am  are  the  corresponding  acceleration  vectors,  we  have 


(4.1) 


(4.2) 


(4.3) 


am  raay  be  measured  by  on-board  missile  sensors.  For  state 
estimation  to  be  used  in  the  missile  guidance  law  development, 
the  target  acceleration  must  be  modeled.  Behavior  of  target 
acceleration  components  is  usually  different  along  target  longi¬ 
tudinal  and  lateral  directions.  To  model  these  components 
separately,  the  targets's  orientation  angle  in  the  inertial  axi 
system  must  be  estimated.  Since  it  is  difficult  to  estimate 
target  orientation  angles,  it  is  desirable  to  model  the  target 
acceleration  components  in  the  same  manner  in  all  directions. 

One  often-used  formulation  is  a  random  walk  model: 


(4.4) 


where  n  is  a  vector  of  white  Gaussian  noise  sources.  Note 


that  the  missile  target  dynamic  model  in  the  inertial  axis  system 
is  linear. 

Nonlinearities  in  the  state  estimator  arise  from  the  measure 
ments.  Passive  seekers  measure  pitch  and  yaw  gimbal  angles  9p 
and  0y  with  respect  to  the  missile  body  frame  (m  refers  to 
quantities  in  the  missile  frame) 

0„  =  arctan  (z  /xm)  ,  (4.5) 

p  mm7 

8  =  arctan  (y  //xf*  +  z^'  )  .  (4.6) 

y  v  m  m  m  i 

The  inertial  line-of-sight  may  also  be  considered  as  potential 
measurements.  The  range  measurement  is 

R  =  (xj  +  yj  +  .  (4.7) 

m  m  m 

Since  the  measurement  model  is  nonlinear,  an  extended  Kalman 
filter  could  be  used  to  estimate  states.  The  covariance  equa¬ 
tions  need  to  be  propagated  in  parallel  with  the  estimated  state 
equations.  Assuming  that  the  missile  velocity  components  may  be 
determined  by  open-loop  integration  of  missile  acceleration 
components,  a  nine-state  Kalman  filter  will  be  required  to  esti¬ 
mate  relative  position,  target  velocity  and  target  acceleration. 
The  covariance  equations  will  place  significant  computation 
requirements  on  the  on-line  processor.  The  time-scale  separation 
approach  will  attempt  to  simplify  this  problem. 

4.2  TIME-SCALE  SEPARATION  FOR  STATE  ESTIMATION 

In  the  state  estimator,  the  time-scale  separation  methods 
will  consider  the  position  as  the  slowest  states,  followed  by 
target  velocity  and  target  acceleration. 

Since  the  measurements  are  all  in  the  slowest  time  frame, 
we  will  attempt  estimation  in  that  time  frame  first.  The  esti¬ 
mator  equations  should  be  of  the  form  (the  superscript  ' " ' 


<  — w  w  w 


denotes  estimate  of  a  quantity): 

X  *  *  * 

X  =  vt  -  vm  +  Kx(y  -  h(x) )  , 
where  the  measurements  are  represented  by 


y  =  h(x)  . 


(4.8) 


(4.9) 


The  gain  K  is  computed  assuming  v+  and  v  doe  not  change 
x  t  m 

significantly  during  this  time  period. 

The  variables  in  the  next  faster  time  scale  will  be  esti¬ 
mated  by  considering  a  measurement  of  vt  based  on  the  position 
equation.  One  possible  model  for  pseudo-velocity  measurement  is 


yv  "  x  +  vm 


■  Kx(y  -  h(x) )  +  vt  . 


(4.10) 


The  estimator  for  target  velocity  components  then  becomes 


t  »  at  +  Kv<yv  '  V 


at  +  Kv  Kx(y  "  h(x))  * 


(4.11) 


In  the  missile  guidance  problem,  noisy  measurements  of  the 
target  position  components  are  available  to  the  missile  in  the 
form  of  look  angles  and  posibly  range.  From  these  measurements 
the  trajectory  of  the  target  is  to  be  estimated.  Further,  this 
trajectory  must  be  updated  with  each  new  measurement.  It  appears 
feasible  to  represent  the  target  acceleration  as  an  auto-regres¬ 
sive  (AR)  model.  Under  the  assumption  that  the  measurement  noise 
has  a  rational  spectrum,  it  can  be  shown  that  the  corresponding 
model  to  be  fitted  to  the  noisy  acceleration  samples,  derived 
from  the  position  measurements,  is  an  auto-regressive  moving- 
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average  (ARMA)  model  whose  denominator  coefficients  are  the 
desired  parameters.  These  parameters  can  be  estimated  using 
recursive  maximum  likelihood  (RML)  or  recursive  prediction  error 
(RPE)  algorithms,  which  are  described  in  Section  4.3.  Once  the 
model  parameters  have  been  estimated,  the  target  position  can  be 
predicted  in  one  of  two  ways:  1)  Estimate  the  acceleration  from 
the  derived  model  and  determine  the  position  estimate  by  inte¬ 
grating  twice,  and  2)  The  target  position  satisfies  an  AR  model 
whose  AR  coefficients  can  be  computed  from  the  above  estimated 
coefficients.  The  target  position  can  then  be  predicted  using 
these  AR  coefficients. 

The  estimated  target  state  variables  serve  as  input  to  the 
algorithm  of  Section  3,  which  computes  the  optimal  missile  state 
variables.  These  nominal  values  are  used  by  an  adaptive  auto¬ 
pilot  as  a  reference;  the  autopilot  provides  actual  control 
signals  to  the  missile  control  surfaces,  in  order  to  track  the 
reference  missile  trajectory  obtained  from  the  optimal  control 
solution.  A  recursive  lattice  algorithm  for  implementing  the 
adaptive  autopilot  is  described  in  Section  4.4. 

4.3  RECURSIVE  MAXIMUM  LIKELIHOOD  (RML)  ARMA  IDENTIFICATION 
Let  the  observed  time  series  be  modeled  as 

A(q_1)  yk  =  C(q_1)  ek 

where 

A(q  )  =  1  a^q  t  ...  4*  Sj  q 

C(q  )  —  1  +  c^q  +  ...  +  Cj^q 

with  q-1  as  the  delay  operator 

q  1  yk  =  yk-l  • 


(4.12) 

(4.13) 
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The  recursive  maximum  likelihood  (RML)  algorithm  computes 
an  estimate  9  of  the  model  parameter  vector  0 


T 

®  (a^  • • •  a^  ^1  • • •  c^) 


by  the  following  set  of  recursions: 


9k  =  9k-l  +  pk  ^k  ek 


where 


pk  ■  pk-x 


pk-i  *11  pk-l 
1  +  *k  *k-i  k 


£k  ^k  -  9i-i  yk 


yk  :=  t"yk-l  *•*  _yk- 


k-L  k-1 


*  ek-NJ 


^k  :=  yk/C(q_1) 


(4.14) 


(4.15) 


(4.16) 


The  Recursive  Prediction  Error  (RPE)  algorithm  uses  modi¬ 
fied  e,  y,  ip  variables.  The  prediction  error  is 

replaced  by  the  a  posteriori  residual 

Ek  :=  yk  "  0k  yk  (4.17) 

where 

yk  :=  t-yk-l  ~yk-L  Gk-1  ek-N^  (4.18) 

The  variable  <l>k  is  also  modified  to 

^k  ^-yk-l  -yk-L’  £k-l  ***  £k-L^  (4.19) 

where 

£k  =  ek/C(q-1)  (4.20) 

A  detailed  discussion  of  both  algorithms  is  provided  in  Appendix 

B. 
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4.4  ADAPTIVE  LATTICE  ALGORITHM  FOR  TRAJECTORY  TRACKING 

The  missile  trajectory  has  to  be  adaptively  controlled  by 
the  autopilot  to  track  the  nominal  trajectory,  which  is  itself 
adaptively  updated  (at  a  lower  rate)  according  to  the  changes 
in  target  trajectory  estimates.  An  algorithm  to  compute  the 
control  inputs  that  achieve  this  objective  is  described  in 
Appendix  D.  The  main  part  of  the  algorithm  (everything  except 
step  V  in  the  Appendix)  is  devoted  to  establish  the  exact  input- 
output  relationship  of  the  missile  and  to  update  this  relation¬ 
ship  as  the  missile  parameters  change.  The  lattice  algorithm 
described  in  Appendix  D  consists  of  a  cascade  connection  of 
identical  lattice  modules  (Figure  4-1).  The  matrix  computations 
performed  by  each  section  can  be  transformed  into  a  set  of  inter¬ 
related  scalar  recursions,  so  that  each  multichannel  module  of 
Figure  4-1  is  replaced  by  a  square  array  of  single-channel  lattice 
cells  (Figure  4-2).  Each  of  these  cells  performs  a  set  of  scalar 
computations  summarized  in  Table  4-1.  The  modular  architecture 
of  the  algorithm  not  only  makes  it  perfectly  suitable  for  VLSI 
implementation  but  also  provides  it  with  a  better  numerical 
behavior  and  higher  throughput  rate  than  any  direct  implementation 
of  the  recursions  of  Appendix  D  on  a  general  purpose  computer. 
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Figure  4-1  Architecture  of  the  Adaptive  Multichannel 
Lattice  Algorithm  for  Trajectory  Tracking 


Figure  4-2  Modular  Architecture  o f  a  Multichannel  Lattice  Section 
e  -  forward  residuals 
r  -  backward  residuals 
,  v  -  auxilliary  signal 


SECTION  5 


SUMMARY 


This  section  summarizes  the  work  performed  under  this 
research  effort.  Plans  for  future  work  are  also  presented. 

5.1  ADVANCED  GUIDANCE  LAWS 

Proportional  navigation  (pronav)  guidance,  which  has  been 
highly  successful  in  the  past,  is  not  likely  to  meet  future 
missile  guidance  requirements;  for  example,  in  an  anti-tactical 
ballistic  missile  (ATBM)  scenario.  A  singular  perturbation 
approximation  to  modern  control  theory  provides  much  high 
accuracy  than  pronav.  The  mathematical  formulation  of  the 
solution  is  significantly  simplified  by  introducing  separation 
of  time  scales.  The  resulting  algorithm  is  well  suited  for 
on-board  real-time  implementation. 

5.2  ADAPTIVE  TARGET  STATE  ESTIMATION  AND  TRACKING 

Efficient  and  fast  algorithms  for  adaptive  target  state 
estimation  and  tracking  have  been  developed.  The  recursive 
maximum  likelihood  (RML)  algorithm  for  fitting  an  auto-regressive 
moving-average  (ARMA)  model  to  noisy  target  position  measurements 
exhibits  fast  convergence  and  is  unbiased  in  the  presence  of 
measurement  noise.  The  adaptive  multichannel  lattice  algorithm 
for  trajectory  tracking  has  excellent  numerical  behavior,  fast 
convergence  and  a  modular  structure  that  makes  it  perfectly 
suitable  for  parallel  processing  implementation. 
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5.3 


FUTURE  RESEARCH 


The  following  areas  need  further  study: 

1.  Real-time  implementation  of  the  advanced  guidance 
law  algorithm  and  evaluation  on  complete  ATBM 
simulation , 

2.  Application  of  the  recursive  maximum  likelihood 
(RML)  algorithm  to  target  trajectory  estimation, 


3.  Application  of  the  recursive  adaptive  lattice-form 
controller  to  design  a  robust  adaptive  autopilot 
for  medium-range  surface-to-air  missiles,  and 

4.  Architectures  for  an  integrated,  parallel,  multi¬ 
microprocessor  implementation  of  the  missile 
guidance  and  control  system. 
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APPENDIX  A 


PROPORTIONAL  NAVIGATION  GUIDANCE 


1.  EQUATIONS 

With  reference  to  Figure  1,  the  engagement  dynamics  in  two 
dimensions  are  [2] 

ro  -  vt  sin  ♦t  -  vm  sin  im  ,  o(tQ)  •  oQ  ,  kD 

r  *  vt  cos  4>t  -  vffl  cos  <*m  ,  r(tQ)  *  rQ  .  (2) 

By  differentiating  the  first  equation,  adding  o  times  the  second 
equation  to  it,  we  get 

ro  +  2ro  +  cos  4>  *  -\r  sin  4  +  cos  <{k0.  +  v  sin  <p.  , 

mm  m  mt  ttt  t 

«<V  *  “o  •  (3) 

Note  that  a  * v  6  .  In  this  section,  we  consider  missiles  with 
m  m  m 

constant  speed,  and  targets  with  no  lateral  acceleration  capabili¬ 
ties  (significant  further  deterioration  occurs  if  these  assumptions 
do  not  hold).  The  change  in  target  speed  is  included  as  follows: 


Figure  1.  Definitions  of  the  Kinematic  Variables 
for  the  Relative  Motion  Between  the  Missile  M 
and  the  Target  T 
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vt  sin  <t>t  . 


(4) 


ro  +  2r o  +  cos  <J>  a  « 

m  m 


The  commanded  acceleration  for  pronav  guidance  law  is  of  the 

form 


amc  "  “An  cos  a  • 
ro 

An  is  the  navigation  gain.  If  the  missile  dynamics  is  negligible, 
i.e.,  am  «  amc,  the  closed-loop  dynamics  is 

ro  +  (2-An)ro  -  vt  sin  $t  , 

which  is  stable  as  long  as  An  >  2  (note  that  r  is  negative). 

Missile  dynamics,  approximated  by  a  second-order  system,  can 
cause  instabilities 


•  2  2 

a.  +  2E  u>  a_  +  u  a_  ■  us  a _ 

m  m  m  m  mm  m  me 


(6) 


A  pronav  guidance  schematic  flowchart  is  shown  in  Figure  2.  The 
characteristic  polynomial  of  the  closed-loop  system  is 


s3* 


/2r  ,  or  \  2  .  /4r  r  ,  2 \ 
\  r  sm  m  f  ^  r  m  m  m J 


+  I 

r 


“m<2' 


V 


(7) 


1.  it. » 


MOMV 

ailOUKE  LAM 
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Figure  2.  Pronav  Guidance  Loop 


2.  ANALYSIS  OF  STABILITY 

The  guidance  law  is  stable,  i.e. ,  the  polynomial  of  Eq.  7 
has  roots  with  negative  real  parts,  if 


Note  that  r  is  nominally  negative  and  -r/r  is  the  time-to-go  if 

there  are  no  maneuvers.  If  A„  exceeds  two,  the  second  and  the 

n 

third  condition  are  always  less  restrictive  than  the  last  condition. 
Therefore,  conditions  for  stability  are 


Note  that  the  pronav  drives  the  missile  unstable  prior  to  impact 

for  all  values  of  damping  ratio.  A  root  locus  plot  with  An  is 

shown  in  Figure  3  to  illustrate  the  stability  problem.  Note  that 

A  >  2  to  stabilize  the  kinematic  pole.  The  point  where  the  complex 
n 

pole  pair  crosses  the  imaginary  axis  depends  on  time-to-go.  As 
the  time-to-go  decreases,  this  crossing  occurs  at  smaller  and 
smaller  A 
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IMAGINARY 

PAR* 


Figure  4  shows  stability  regions  for  various  pronav  gains 
and  two  values  of  missile  damping.  The  higher  the  pronav  gain, 
the  earlier  the  missile  goes  unstable. 

3.  SUMMARY 

When  the  target  has  high  acceleration  components  the  terminal 
instability  of  pronav  causes  large  miss  distances.  The  pronav  is 
also  too  slow  to  respond  to  large  errors  when  the  range  is  large, 
causing  missile  saturation  prior  to  impact. 


4 

PRONAV 

CAIN 


5C 


TIW.TO.QO  (UNITS  OP  UmJ 

Figure  4.  Regions  of  Missile  Stability 


APPENDIX  B 


A  QUASI -NEWTON  ALGORITHM  FOR  ML  ESTIMATION 
OF  ARMA  MODEL  PARAMETERS 


Let  the  observed  time  series  (y(t),  t  *  0,  1, 
modelled  as 


where 


. . }  be 


A(q-1) 

y(t)  ®  c(q_1)  e(t ) 

(1) 

ACq"1) 

*  1  +  a^q-1  +  . . .  +  aLq~L 

C(q-1) 

=  1  +  c^q-^  +  .  .  .  +  Cj^q-^ 

(2) 

with  q  1  as  the  delay  operator 


-1 


q  y(t)  =  yCt-1)  , 

and  e(t )  is  a  sequence  of  zero  mean  white  noise  samples. 
Rewriting  (1)  as 


y(t)  = 


1  .  ACq'1) 
C(q‘l) 


y(t)  +  e(t) 


the  one-step  ahead  predictor  y < 1 1 6 )  is  given  by 
y(t | e)  = 


i  -  m<l±> 


C(q-1 ) 


y(t ) 


(3) 


(4) 


where  ■:  is  the  parameter  vector  as  defined  in  (6), 
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From  (3)  and  (4),  the  prediction  error  e(t,5)  is 


e(t,5)  =  y (t )  -  y(t|6)  =  y(t)6T;(t,6) 


where 


@  *  ( a  ^  •••  aL  c ^  «  •  •  c )  ( 6 ) 

CT(t,e)  »  (-y(t-l)  ...  -y(t-L)  eCt-l,e)  ...  e(t-N,9)) 

The  off-line  maximum  likelihood  (ML)  method  of  estimating 
the  parameter  vector  6  corresponds  to  minimizing  the  function 


J  -s 


E*(s,e> 


The  above  minimization  problem  is  non-linear  in  9  and  hence 
an  explicit  solution  is  not  possible.  Therefore,  a  numerical 
search  procedure  based  on  a  Quasi-Newton  method  will  be  used  to 
find  9. 

Differentiating  J  with  respect  to  0  gives 


VJ  « 


e (s , 6)  Ve(s,0) 


where  Vj  and  Ve(s,9)  denote  the  gradient  vector  of  J  and 
e(s,a)»  respectively,  with  respect  to  9.  Differentiating  once 
again  gives 


[ Vc (s , 0 ) VTe (s , 0 )  +  V2e(s,e)e(s,0)] 


From  the  expression  (5), 


Ve(s,9)  *  -  Vy(s|0)  ■  -  ^(s,0) 


(10) 
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where  ^(s ,r>  denotes  the  gradient  vector  of  the  prediction 
y  - '  r  ) . 

2 

The  Hessian,  7  J,  will  now  be  approximated.  At  the  true 
minimum  of  (8),  the  second  term  of  (9)  can  be  shown  to  be  zero. 
Since  the  true  Hessian  is  more  important  close  to  the  true  mini¬ 
mum  than  elsewhere,  we  approximate  the  Hessian  by  the  first  term 
of  (9)  which,  in  view  of  (10),  becomes 


V2J 


'V 


^(8,6)  (s, 9) 


(11) 


Further  approximations  are  needed  in  order  to  obtain  a  recursive 
Quasi-Newton  algorithm  from  (8)  and  (11).  Computation  of  e(t,9) 
requires  all  the  data  up  to  t.  This  computation  is  approximated 
by  using  the  latest  values  of  the  data  and  the  parameter  estimates, 
and  denote  the  correponding  C(t,8),  iKt,8)  and  e(t,8)  by 
C(t),  4>(t)  and  e(t),  respectively.  The  Quasi-Newton  update 

of  the  parameter  vector  is  then  given  by: 


9  (t ) 


8(t-l) 


i-l 


^  Cs ) 


*T(s) 


4Kt)  e(t) 


(12) 


where  we  approximated  VJ  by  ty(t)e(t).  Note  that  the  effect 
of  the  above  approximations  on  the  asymptotic  values  of  the 
parameter  estimates  is  negligible  if  the  roots  of  C(z)  lie 
inside  the  unit  circle. 


Now  consider  the  computation  of  the  gradient  vector  of  the 
prediction.  From  (4), 


y(t ) 

which  gives 


!  _  Aialil 

CCq'1) 


y(t ) 


y(t-i) 

C(q_1) 


a, 

-  y(t-i) 


(13) 
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and 


3y(t) 

3ci 


e(t-i) 
C (q-1 ) 


c(t-i) 


Combining  (3)  and  (14),  we  obtain 


(14) 


Vy(t)  =  ^ (t )  =  .  (15) 

C (q-1) 

<v 

y(t)  and  e(t)  are  the  filtered  variables.  The  vector  ^(t) 
can  be  defined  in  terms  of  these  variables  as 


^T(t)  *  [-y(t-l)  ...  -y(t-L)  E(t-l)  ...  E  ( t -N ) ]  (16) 

In  computing  4'(t),  the  estimated  value  of  6  at  t-1,  6(t-l), 
is  used.  It  is  then  easy  to  see  that 


'V 

y(t ) 

=  y(t ) 

-  Cj(t )y(t-l)  . . . 

-  cN(t) 

y(t-N) 

(17) 

€(t) 

=  c(t) 

-  cx(t )e(t-l)  . . . 

-  cN(t) 

E(t-N) 

(18) 

Using  the  matrix-inversion  lemma,  a  recursive  version  for 
(2)  can  be  obtained.  The  resulting  algorithm  follows: 


A  As 


6(t ) 

*  e (t-i ) 

+  P(t)  4(t)  e(t ) 

P(t) 

*  P(t-i) 

P(t-l)  *(t)  \l?( t) 

P(t-l) 

1  +  *T(t)  P(t-l) 

'l(t) 

E  ( t  ) 

=  y (t )  - 

eVt-l)  ;(t) 

01  (t) 

=  (a1(t) 

•  •  •  (t  )  Cj(t  )  •  • 

A 

•  Cjq  ( t  ) 

C(t)  and  ij<(t)  are  as  defined  in  (6)  and  (16),  respectively.  This 
algorithm  is  called  the  recursive  maximum  likelihood  (RML)  method. 
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Note  here  that  the  latest  error  term  contained  in  £(t)  is 
i ( t - 1 ) ,  and  the  latest  filtered  error  in  ;(t)  is  e(t-l). 

A 

Since  the  estimate  6(t-l),  available  at  the  beginning  of  the 
t-th  sampling  interval,  facilitates  the  computation  of  the  a 
posteriori  residual  ?(t-l),  defined  as 

e(t-l)  =  y(t-l)  -  9T(t-l)  s(t-l>  ,  '  (20) 

the  residuals  can  be  used  in  place  of  the  prediction  errors  in 
C(t)  and  i(t) .  Thus,  the  modified  forms  of  ?(t)  and  iHt) 
are  given  by 

CT(t)  *  (-y(t-l)  ...  -y(t-L)  F(t-l)  ...  e(t-N)) 

(21) 

VT(t)  *  (-y(t-l)  ...  -y(t-L)  e(t-l)  ...  e(t-N)) 

—  —  -1 

where  e(t)  =  e(t)/C(q  ).  The  algorithm  (19)  with  ;(t)  and 
_,(t)  as  defined  in  (21)  is  called  the  recursive  prediction  error 
method  ( RPEM) .  The  difference  between  (6),  (16),  and  (21)  re¬ 
flects  in  the  transient  behavior  of  the  RML  method  and  the  RPEM; 
however,  their  asymptotic  behavior  is  identical. 
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APPENDIX  C 


THE  ADAPTIVE  LEAST-SQUARES  PROBLEM 

Exact-least-squares  lattice  algorithms  provide  recursive 
solutions  to  the  following  adaptive  least  squares  problem: 

Given  two  sequences  of  multichannel  measurements 

y(0) ,  y(l) , . . . 

and 

x(0) ,  x(l) , . . . 

find  a  linear  estimate  of  x(k)  based  on  m  previous 
measurements  of  y  ,  namely 
m 

x(k)  :*  £  h.y(k-i) 
i=l  1 

such  that  the  exponentially  weighted  cost  function 

C®’*  :=  E  Xt_k  |  |x(k)-x(k)|  |2 
T  k=0 

is  minimized. 

The  optimal  solution  to  this  problem,  is  a  function  of  m,X,  t 
and  also  of  the  data  y(k),  x(k).  Since  it  is  customary  in 
adaptive  applications  to  solve  the  problem  for  several  values  of 
ro  and  t  ,  the  dependence  of  the  solution  upon  these  parameters 
will  be  indicated  explicitly  by  introducing  the  notation 


for  the  solution  that  minimizes  Cm'*  .  The  exponential  weight 
is  usually  chosen  in  the  range  0.9  <  1  <  1.0. 
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P 

f. 

[< 


I 
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Adaptive  lattice  algorithms  provide  solutions  which  are 
recursive  both  in  the  parameter  t  ("time”)  and  in  the  parameter 
m  ("order").  These  algorithms  are  based,  consequently,  on  two 
update  formulas:  an  order-update,  which  relates  the  solution 

h1?-?"  to  h1!1  .  ,  and  a  time-update,  which  relates  h"1  1  to 

h™  .  .  These  update  formulas  can  be  expressed  either  in  terms  of 

i » t 


the  residuals 


e  (t)  : 
m 


=  x(t) 


m 

-  E 

i=  1 


hi,t 


y(t-i) 


(l) 


or  in  terms  of  the  predict ion -errors 

ES(t>  x<t>  -  i/T.t-i 


(2) 


The  derivation  of  the  time-  and  order-update  formulas  is  based 
upon  a  geometric  approach  which  is  described  in  the  sequel. 


Define  the  rectangular  matrices 


yt  :=  [y(t),  y(t-i) , . . . ,  y(i),  y(0),  0...0] 
xt  :«=  [x(t),  x(t-l),...,  x(l),  x(  0) ,  0...0] 


(3) 


The  number  of  rows  in  each  matrix  is  determined  by  the  number  of 
channels  of  the  corresponding  signal.  The  length  of  each  row 
(say  N  )  is  chosen  large  enough  to  guarantee  the  appearance  of 
zero  columns  of  the  end  of  x^ ,  y  for  every  t  in  consideration. 
The  rows  of  x^.y.j.  are  ‘t^ere:^ore  elements  (vectors)  in  the  linear 
space  of  row  vectors  of  length  N  .  Defining  an  inner  product  for 
every  two  row  vectors  of  length  N, 


<a,b>  :=  aAb* 

A  :*  diag  { >  * ;  0  £  i  <_  N-l) 


(4) 


where  the  asterisk  *  denotes  the  Hermitian  transpose,  we  obtain 
a  (weighted)  Euclidean  space.  A  collection  of  several  vectors  in 
this  space,  say  x1,x2,...,xk  ,  written  formally  as  a  column 
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I 


*  -(=0 


(5) 


will  be  called  a  vector  array  (VA). 
VAs  X,Y  is  defined  as  the  matrix 

<X,Y>  :=  [<xi,yi>] 


The  (inner)  product  of  two 


(6) 


whose  (i,j)  element  is  <xi,yi>,  the  inner  product  of  the 
i-th  vector  in  X  with  the  j-th  vector  in  Y  . 

The  cost  function 


can  be  expressed  in  terms  of  the 


VAs 


xt’yt 


as 


Cm> t  =  tr  <  xt 


m 


m  m 

-  £  Vt-i  •  xt  -  £  Vt-i’ 

i=l  i=l 

that  minimizes  the  cost  function  is  obtained, 
(i.e.,  projecting  every  vector  in  x.  ) 


The  solution  Ik  t 

by  projecting  the  VA  xt  - -  - „ - ,  - -  „ 

on  the  subspace  spanned  by  all  the  vectors  contained  in  the  VAs 
yt_t,  yt_m  .  Since  projections  play  a  central  role  in 

solving  the  least  squares  problem,  it  will  be  convenient  to  in¬ 
troduce  a  shorthand  notation.  Let  x  denote  the  residual  ob¬ 
tained  by  subtracting  from  every  vector  xi  in  the  VA  X  ,  the 
projection  of  xA  on  the  subspace  spanned  by  the  vectors  con¬ 
tained  in  U 
as 


Then  the  minimal  cost  can  therefore  be  expressed 


min  C 


,m,  x  _ 


m 


t  =  tr  <em,t ’  Em,t> 

.  : *  ( x.  ) ,  .  =  x,  -  h™ 

m't  1  (yt-l . yt-m>  1  ~ 


(7) 


t  yt-i 


where  e  *  is  the  residual  of  x*.  after  removing  its  projec- 
m ,  t  t 

tion  on  span  {y.  ,  y*  _}  . 

The  order-  and  time-update  formulas  of  exact-least-squares 
lattice  algorithms  involve  only  the  first  column  of  the  matrix 
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,  t  ;  namely. 


where 


t )  =  <e  +  ,  ir> 
m  m ,  t 


TT  :  =  [  1,  0,  . .  . ,  0  ] 


Notice  that  em(t)  is  precisely  the  multichannel  residual  defined 
by  (1).  The  prediction-error  ej^(t)  »  defined  by  (2),  is, 
similarly,  the  first  column  of  the  matrix 


i  t  :=  xt  "  1]  hi  t-1  yt  i 

’  i=l 


(10) 


which  seems  to  have  no  geometric  interpretation  in  terms  of 
projections.  If  the  columns  of  this  matrix  are  shifted  one  step 
to  the  left,  and  a  column  of  zeros  is  introduced  at  the  extreme 
right,  the  resulting  matrix  is 

m 

Xt-1  *  £  hi,t-l  yt-l-i  5  em, t-1 

i«l 

which  is,  of  course,  the  residual  of  projecting  xt  on  span 
(y+  o  y+  „  - }  .  Denoting  the  left-shift  on  row  vectors  by 

L  f  •  •  •  f  x  "in*  j. 

D  we  have 


Dyt  -  y, 


(11) 


and  as  a  consequence  of  our  last  argument 

e^(t)  *  <eE  +  »  v> 

ID  ID,  t 

*  c  ^ 
m, t  m, t-1 


(12) 


This  identity  can  be  used  to  derive  order-  and  time-update  for¬ 
mulas  in  terms  of  prediction-errors  rather  than  residuals. 
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APPENDIX  D 

ADAPTIVE  CONTROL  AND  PREDICTION  USING  LATTICE  STRUCTURES 

\ 

A  large  number  of  adaptive  control  and  prediction  algorithms 
can  be  described  exactly  as  a  lattice  form  implementation  that 

i.  u 

requires  only  0(M)  computations  per  time  update  for  an  M  order 
model.  The  well-known  advantages  of  lattice  form  implementation 
in  terms  of  numerical  stability  and  simultaneous  availability  of 
lower  order  models  become  available  for  algorithms  with  known, 
desirable  asymptotic  convergence  properties. 

For  the  ARMAX  (Autoregressive  Moving  Average  with  Exogenous 
Inputs)  models,  current  state-of-the-art  adaptive  prediction 
either  risks  convergence  to  a  local  maxima  of  the  likelihood 
function,  or  requires  a  Strict  Positive  Real  condition  for  con¬ 
vergence.  Over-parameterization  of  the  predictor  can  be  used 
to  ensure  that  neither  of  these  two  problems  would  arise  nor 
would  the  convergence  rate  reduce  substantially.  This  is  attrac¬ 
tive  only  because  of  the  0(M)  computational  complexity  of  the 
algorithm. 

INTRODUCTION 

Consider  a  linear  system  with  p  inputs,  p  outputs  described 

by: 

A(q_1)yt  =  q"1B(q_1)ut  +  C(q_1)vt  (I) 

q-1  =  delay  operator  in  discrete-time.  {yt>,  {ut),  and  (vt) 

denote,  respectively,  the  output,  the  input  and  the  noise  process, 
each  of  dimension  pxl . 

ACq'1)  -  I  ♦  Ajq'1  ♦  ...  +  Anq-n 
C(q_1)  =  I  +  Cjq"1  +  ...  +  cnq’n 
B(q-1)  *  Bq  +  Bjq-1  +  ...  +  Bmq"m 
det  (Bq)  #  0  . 
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Objective  of  the  adaptive  controller  is  to  make  the  output  follow 
a  reference  signal  yreft  in  minimum_variance  sense.  Viz. 

MiD  E(|yt+1  -  yreft+il2|yt.ut.yreft+i,yt.r»t.1,yreft,...} 

Equation  (I)  can  be  rewritten  in  its  prediction  error  form  as 

C(q_1H(yt+1  -  yreft+1>  '  vt+l}  “  4 

4  (C(«'1)  -  AO-1))  •  <yt  -  yreft>  -  A(Q-1)yrej  } 


*  Vt  -  W 


where 


=  tut’ut-l»  *  **’ (yt"treft)T,(yt-l“yreftl)  ’  **•’ 

T  T  -.T 

yreft’  yreft_^’  ** ' 

Let  r  =  max( (n-1) ,m) ,  dim  <4t  *  3*(r+l)4pxl  and  6Q  is  a  matrix 
of  coefficients  with  dimension  *  p  x  dim  .  0Q  may  have  zero 
elements.  Model  order  M  =  r  +  1. 

6q  —  [  Bq  ,  , B2 1  •  •  « » By  ,  (Cj —A.^ )  >  ( ^2~^2  ^  *  *  *  *  *  (Cr+i ^  * 

-  Aj,  -  Ag »  ...»  -  Ar+1 ] 

Recursive  least  squares  algorithm  for  parameter  update  in  the 
regression  form 


(yt*l  -  yref  >  4  yref. 


-  6*t  . 


(ID) 


is  used  to  give  parameter  estimate  w 

The  parameter  estimate  is  used  in  turn  to  compute  the  control 
ut  in  the  following  manner: 


Choose  ut  such  that  ©t  -  yref  *  0 


(CONTROL) 
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Such  a  choice  of  ut  gives  reference  following  in  the  minimum- 
variance  sense. 


ID  equations  are  usually  implemented  as: 


-  e. 


+  p. 


<x+  <t>: 


- 1  ~  t-i  t-i  t-iv~"  vt-i1  t-iyt-i '  yjrt  -  “t-i^t-i' 


P 


t 


P 


t-1 


Pt-l^t-l^t-lPt-l 
x  +  *t-lPt-l*t-l 


1 

X 


Control  u^  is  computed  by  solving  the  linear  set  of  equations 


t*t 


ref 


t+1 


Note  that  this  involves  inversion  of  a  pxp  matrix  consisting  of 
the  first  pxp  matrix  of  0  . 

This  algorithm  of  (1)  finding  a  linear  least  squares  esti¬ 
mate  in  a  linear  regression  model  and  (2)  then  computing  the 
control  to  make  the  prediction  based  on  the  current  parameter 
estimates  equal  to  some  known  value,  lies  at  the  heart  of  a  large 
number  of  "successful"  adaptive  control  algorithms.  The  same 
basic  algorithm  is  also  used  for  recursive  prediction  algorithms. 


Lattice  Form  Algorithm 


Exact  implementation  of  a  recursive  least-squares  algorithm 
can  be  done  using  the  joint-process  ladder  form  with  pre-window 
ing.  The  parameter  estimate  is  In  terms  of  the  reflection  co- 

A 

efficients  and  not  in  terms  of  0.  The  reflection  coefficients 
are  used  to  compute  the  prediction  and  efficiently  obtain  the 
control  ut  to  make  the  predicted  value  equal  to  the  desired 
value . 


The  proposed  procedure  involves  p+1  iterations  of  computing 
prediction  ol  <yt+1-yrelt+i>  for  ut  -  0,  ,  . . p.  «t 
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1 


denotes  the  ith  Euclidean  basis  vector  [00.  . 
Each  such  prediction  involves  0(M)  operations, 
linear  affine  equations  in  u.  is  then  solved  in 


.  .  0]  in 
The  set  of 

0(p3> 


EP. 


operations 

to  obtain  ut .  The  total  operation  count  thus  remains  0(M),  i.e., 
linear  in  the  model  order.  Note  that  since  we  are  using  the 
ladder  form  merely  to  implement  the  exact  equivalent  of  the  algo- 

A 

rithm  using  explicit  estimation  of  6,  the  convergence  properties 
and  stability  property  of  the  original  algorithm  continues  to 
hold  for  the  ladder  case.  Let  us  describe  the  algorithm  in  detail 
below: 


(I)  Initialize  at  time  t  *  0;  n  *  0,  ...,  M-l 

R~e(0)  =  61 
R~r(-1)  «  61 
Fn(0)  =  0 
rn(0)  *  0 

yn<0)  *  o 

(II)  At  time  t  we  have  in  memory  [R”e(t-l)]t  F  (t-1), 

Fj(t-l),  yn<t>,  [R^r(t-2)],  rn(t-l),  Kj'(t-l) 

Compute:  n  =  0,  ...,  M-l 

♦(Remark  on  notation:  The  subscript  n  denotes  the  lattice 
stage.  To  avoid  confusion,  the  time  parameter  t  is  now 
appearing  in  (  ) . ) 


Kn(t-1)  *=  F^(t-l)[R”r(t-2)] 
K*(t-1)  =  Fn(t-1) [R~e(t-1)] 


R~r(t-1) 


„  R”r(t-2 )  r  (t-l)r^(t-l)R”r(t-2 ) 

R~rf*  n _ n _ n _ n _ 

11  \  ♦  A  /^  \  N  .-,!/*■  i  Vt5  1 


X(l-Sn(t))+r*(t-l)Rni(t-2)rn(t-l) 


Bn+l(t)  *  6n(t)  +  <l-^nCt))2r^(t-l)R^r(t-l)rn(t-l) 


1 

X 


e0(t)  -  o 
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(Ill)  Measurement  y(t)  is  used 

en(t)  =  y(t )  -  yn(t)  n  *  0,  M-l 

Fn(t)  "  ^t-1)  +  d-6n(t))rn(t-l)ej(t) 

Kn(t)  =  Fn(t)tRnr(t"1)] 


(IV)  Let  *u(t) 


ei 


i  ~  1  *  •••»  P 


Let  °u(t)  =  0 

Let  ix(t)  =  [1u(t)T,  yT(t),  y^ef(t+l)]T 
Compute  for  i  =  1,  p;  n  =  1»  ...»  M 

*en(,)  *  len-l(t>  + 

V*’  *  + 

t+n  -  -  uLifo'viC) 

ieQ(t)  =  lr0(t)  *  Xx(t) 
giving  p  predictions  yM(t+l) 

(V)  Solve  the  set  of  linear  equations  defined  by 


where 


Au(t)  =  b 

»  -  yre,<t+n  -  Vt+1> 


and 

column  i  of  A  =  (*yn(t+l)  -  b) 

This  choice  of  u(t)  will  give 

yM(t+l)  ■  yref (t+1)  as  desired. 

(VI)  Apply  control  u(t). 
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(VII) 


Use  x(t)  formed  by  u(t),  y(t),  yref(t)  and  compute 

for  n  -  l,  m 

eQ(t)  -  rQ(t )  -  x(t ) 

en(t>  *  ♦  Vl<t-1>rn-l<t-1) 

VtJ  -  +  ClO-'l'n-l1’) 

yn<t+l>  *  W"11  -  Kn-l(t>rn-l(t:> 


» (l-6n(t )  )+e*(t  )B"e(t-l)en(t ) 


T,  ^  lri-e, 


B;e(t-1)  - 


Fn(t)  -  XFn(t-l)  +  (l-8n(t))rn(t-l)ej(t) 
Go  to  step  II 
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